#!/bin/bash
PANF=/home/cmbpanfs-01/yanghoch # can be overwritten by bashrc
source ~/.bashrc
source ~/SHELL/isFileEmpty.sh
if [ $# -ne 4 ]; then
	echo "This shell enumerate isoforms and estimate their expressed levels"
	echo "in a given gene region. It uses:"
	echo "(1) the sequence (.fa)"
	echo "(2) the low coverage region"
	echo "(3) the junctions (.sam + .junctions ->.bed)"	
	echo "It also generates the coverage .bed file of the given region"
	echo ""
	echo "The 1st ptr is gene name"
	echo "The 2nd ptr is chr seq name"
	echo "The 3rd ptr is the range start"
	echo "The 4th ptr is the range end"
	echo "4 ptrs are needed"
	echo "$# is provided" 
	exit 0
fi

# use gene name as working directory
if ! test -d $1; then
	mkdir $1
fi
cd $1
echo "Working in the $PWD"
# create tmp directory 
TMP_DIR='./tmp'
if ! test -d $TMP_DIR; then
	mkdir $TMP_DIR
fi

DATA_DIR=$PANF/'hg18'
SEQ_PATH=$DATA_DIR/"$2.fa"        # The whole sequence
SEQ_OUT=$TMP_DIR/"$1.fa"          # The extracted sequence
SEA_PATH=$DATA_DIR/"$2.sea"       # The whole low coverage regions 
SEA_OUT=$TMP_DIR/"$1.sea"         # The extracted low coverage regions
SAM_PATH=$DATA_DIR/"$2.sam"       # The whole junctions support mappings
SAM_OUT=$TMP_DIR/"$1.sam"         # The extracted detected junctions
SAM2BED_JUC=$TMP_DIR/"$1.bed"     # The extracted knwon junctions
JUC_PATH=$DATA_DIR/"$2.junctions" # The whole known junctions
JUC_OUT=$TMP_DIR/"$1.junctions"   # The extracted known junctions  

READS_DIR=$PANF/'BrainTran/'
#READS_PATH=$READS_DIR/'brain_seq.fastq'
#READS_PATH=$READS_DIR/'brain_seq.txt'

# The path for the list with (1)IM_READS, (2)IM_JUNC_READS and (3)IUM_READS
READS_PATH=$TMP_DIR/$1.txt
# The path for the initially mapped reads
IM_READS_PATH=$TMP_DIR/'MapBrainSeq.fa' 
# The path for the initially mapped result, sam format -> IM_READS.fa
IM_SAM_PATH=$READS_DIR/'Mapping/'$2.map
# The path for the initially junction support reads
IM_JUNC_READS_PATH=$TMP_DIR/'MapJuncBrainSeq.fa'
# The path for the junction support mapping result, sam format -> IM_SAM_JUNC_PATH  
IM_SAM_JUNC_PATH=$SAM_OUT
# The path for the initially unmapped reads
IUM_READS_PATH=$READS_DIR/'unMapBrainSeq.fa' 

# extract the initially mapped reads
python ~/SCRIPT/extract_mapped_reads_from_sam.py\
 $IM_SAM_PATH $IM_READS_PATH $3 $4
# extract the initially mapped junction-support reads
python ~/SCRIPT/extract_mapped_reads_from_sam.py \
 $IM_SAM_JUNC_PATH $IM_JUNC_READS_PATH 0 $(($4 - $3 + 1))

#COVERAGE_THRESHOLD=1
#NO_MAPPED_READS=$(wc -l $IM_READS_PATH | cut -f 1 -d' ')
#if [ $NO_MAPPED_READS -le $(($COVERAGE_THRESHOLD * $(($4 - $3 )))) ]; then
#	\rm $IM_READS_PATH
#	echo "$1 has too low coverage" 
#	exit
#fi
# create the list of the IM and IUM reads 
echo "$IUM_READS_PATH" > $READS_PATH
echo "$IM_READS_PATH" >> $READS_PATH
echo "$IM_JUNC_READS_PATH" >> $READS_PATH

a=$(isFileEmpty "$SEQ_PATH")
b=$(isFileEmpty "$SEA_PATH")
c=$(isFileEmpty "$SAM_PATH")
d=$(isFileEmpty "$WIG_PATH")
e=$(isFileEmpty "$READS_PATH")

echo ""
echo "\1. Extract seq in $SEQ_PATH from $3 to $4."
python ~/SCRIPT/getsequences.py $SEQ_PATH $SEQ_OUT $3 $4

echo ""
echo "2. Extract low-coverage regions and junctions in the gene range."
SHIFT=$3
# Get low coverage from (.sea) .wig file
python ~/SCRIPT/one_contig_one_oceans_file.py\
 $SEA_PATH $SEA_OUT $3 $4 $SHIFT
# Get junctions from sam file
python ~/SCRIPT/one_contig_one_junctions_file.py\
 $SAM_PATH $SAM_OUT $3 $4 $((SHIFT+1)) 
# Get known jnctions from bed file
BED_RANGE_COL=1
python ~/SCRIPT/search_a_large_sorted_file.py\
	 $JUC_PATH $JUC_OUT $BED_RANGE_COL $3 $4 $SHIFT
# change sam file to bed file for junctions,
python ~/SCRIPT/one_contig_one_junctions_file.py $SAM_OUT $SAM2BED_JUC $2
if ! test -s $SAM2BED_JUC ; then
	echo "No detected junction!"
	exit
elif ! test -s $JUC_OUT ; then
	echo "No known junction!" 
else		
	### currently marked out because only Tophat junctions are used.
	# merge with known junctions to BED
	MERGE_TMP=$1.merge
	~/SHELL/mergeBED.sh $SAM2BED_JUC $JUC_OUT $MERGE_TMP	
	#mv $MERGE_TMP $SAM2BED_JUC

	# detect new junctions
	# Substract the known junctions from detect junctions to find new junctions
	~/SHELL/except.sh 1_$(basename $SAM2BED_JUC) 2_$(basename $JUC_OUT)
	#~/SHELL/except.sh $SAM2BED_JUC $JUC_OUT
	# Create pattern files to see the sequence of new junctions.	
fi 

echo ""
echo "3. Get the isoform, CDS, mapped reads, and the abundance levels"
GENE=$TMP_DIR/"$1"
~/SHELL/PerM_TranSEQ.sh $GENE $READS_PATH
echo ""
echo "4. Get the html for BED coverage and the most abundant isoforms"
if test -s $GENE'_I.fa'; then 
	~/SHELL/outRNAseqHTML.sh $1 $2 $3 $4
else 
	echo "No valid isoform; Skip from HTML report"
fi
## clean up
exit #tmperariliy don't clean up
\rm *.log
echo "tmp dir $TMP_DIR"
if [ $TMP_DIR == "./tmp" ]; then
	echo "Delete tmp dir"
	for file in $(ls ./tmp); do
		\rm ./tmp/$file
	done
	\rmdir "./tmp"
fi
exit
